function h = spektrHVLn(q,n,atomic_number)

%%**************************************************************************
%% System name:      SPEKTR
%% Module name:      spektrHVLn.m
%% Version number:   1
%% Revision number:  00
%% Revision date:    15-Mar-2004
%%
%% 2016 (C) Copyright by Jeffrey H. Siewerdsen.
%%          I-STAR Lab
%%          Johns Hopkins University
%%
%%  Usage:  h = spektrHVLn(q, n, atomic_number)
%%
%%      eg. h = spektrHVLn(q, 1, 13)
%%      calculates the first HVL of Aluminum
%%
%%  Input Parameters:
%%      q - x-ray energy spectrum (1-150keV)
%%      n - number of the half value layers to calculate
%%      atomic_number - element comprising the filter of interest.
%%
%%  Output Parameters:
%%      h - thickness of n'th HVL (Half Value Layer) [mm]
%%
%%  Description:
%%      This function calculates the n'th HVL (Half Value Layer),
%%      which is the thickness required to attenuate the beam
%%      intensity by 1/2^n.
%%
%%      This function accounts for the polyenergetic case.
%%
%%  Notes:
%%
%%*************************************************************************
%% References:
%%
%%*************************************************************************
%% Revision History
%%  0.000    2003 05 01     AW  Initial code
%%	1.000    2004 03 15     DJM Initial released version
%%*************************************************************************
%%

%%% ... accuracy of search algorithm
accuracy = 0.0001;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial thickness of element
thickness = 0;
min = 0;
max = 100;

% calculate the thickness in which the exposure to due varying spectrum is
% reduced to 1/2^n  of it's initial exposure

exp_temp = spektrExposure(q);
exp_target = (2^-n) * exp_temp;

q_min = spektrBeers(q,[atomic_number, max]);
exp_min = spektrExposure(q_min);
while exp_min > exp_target
    min = max;
    max = 2 * max;
    q_min = spektrBeers(q,[atomic_number, max]);
    exp_min = spektrExposure(q_min);
end

while(exp_temp > exp_target * (1 + accuracy)) || (exp_temp < exp_target * (1 - accuracy))

    %incrementing the thickness of the filter

    thickness = (min + max)/2;

    %   generate spectrum given a elemental filter of a specified thickness
    q_temp = spektrBeers(q,[atomic_number, thickness]);

    %   calculate the exposure of the spectrum
    exp_temp = spektrExposure(q_temp);

    if exp_temp < exp_target
        max = thickness;
    elseif exp_temp > exp_target
        min = thickness;
    end

end

% Return HVL
h = thickness;